/* 2010: Copyright by Massimo Bernaschi, Giorgio Parisi, Leonardo Parisi */
#include <time.h>
#include <cutil.h>

#include "hsggpuSM4P.h"
#include "timing.h"

#include "dictionary.h"
#include "iniparser.h"
#include "util.h"

#include "MersenneTwister.h"
#if !defined(HOSTREDU)
#include "threadFenceReduction_kernel.cu"
#endif


//-------------------------------------------------------------------------------------------------------
// Usage
//-------------------------------------------------------------------------------------------------------
void Usage(char *program) {
  fprintf(stderr,"Usage: %s -i inputfile [-r spin_file] [-o (over-relax only)] [-b (heat bath only)]\n",program);
}


//-------------------------------------------------------------------------------------------------------
// Init
//-------------------------------------------------------------------------------------------------------
void Init(void){

  if(sizeof(REAL)==sizeof(float)) {
    zero=0.0f;
    one=1.0f;
    two=2.0f;
    three=3.0f;
    zerofive=0.5f;
    zerotwentyfive=0.25f;
  } else if(sizeof(REAL)==sizeof(double)) {
    zero=0.0;
    one=1.0;
    two=2.0;
    three=3.0;
    zerofive=0.5;
    zerotwentyfive=0.25;
  } else {
    writelog(TRUE,APPLICATION_RC,"Invalid size of REAL: %d",sizeof(REAL));
  }

  zeroEne = 0.0; //MESSA SOLO PER IL CALCOLO DELL'ENERGIA IN DOUBLE COME PROVA

  h_a = (REAL*) Malloc(ms_s);
  h_b = (REAL*) Malloc(ms_s);
  h_c = (REAL*) Malloc(ms_s);

  h_Jpx = (REAL*) Malloc(ms_j);
  h_Jpy = (REAL*) Malloc(ms_j);
  h_Jpz = (REAL*) Malloc(ms_j);

  h_Energy = (REAL*) Malloc(ms_s);

  h_sdbL = findlog2(h_side); //log_2(L) used to divide by L
  h_sdbLs = findlog2(h_side*h_side); //log_2(L*L) used to divide by (L*L)

  h_maskL = h_side-1;
  h_maskL2 = (h_side*h_side)-1;
  h_maskL3 = V(h_side)-1;
  h_Ldb2 = h_side>>1;
  h_Lsdb2 = (h_side*h_side)>>1;

  h_nspin=V(h_side);

  srand48(seed);

}


//-------------------------------------------------------------------------------------------------------
// CopyConstOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyConstOnDevice(void){

  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_side",&h_side,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_nspin",&h_nspin,sizeof(int),0,cudaMemcpyHostToDevice) );

  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_sdbL",&h_sdbL,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_sdbLs",&h_sdbLs,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_maskL",&h_maskL,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_maskL2",&h_maskL2,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_maskL3",&h_maskL3,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_Ldb2",&h_Ldb2,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_Lsdb2",&h_Lsdb2,sizeof(int),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_h",&h_h,sizeof(REAL),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_Beta",&h_Beta,sizeof(REAL),0,cudaMemcpyHostToDevice) );

  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_const_one",&one,sizeof(REAL),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_const_two",&two,sizeof(REAL),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_const_zerotwentyfive",&zerotwentyfive,sizeof(REAL),0,cudaMemcpyHostToDevice) );

  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_overrelaxhits",&overrelaxhits,sizeof(int),0,cudaMemcpyHostToDevice) );

  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_const_shmpervar",&shmpervar,sizeof(shmpervar),0,cudaMemcpyHostToDevice) );
  CUDA_SAFE_CALL( cudaMemcpyToSymbol("d_const_shmnitemspp",&shmnitemspp,sizeof(shmnitemspp),0,cudaMemcpyHostToDevice) );

}


//-------------------------------------------------------------------------------------------------------
// findlog2
//-------------------------------------------------------------------------------------------------------
int findlog2(int n){

  int i=0;

  if(n<1 || n&1) { writelog(TRUE,FINDLOG_RC,"Invalid value %d. Must be > 0 and even\n",n); }

  while(n>0) {
    n>>=1;
    if(n>0) i++;
  }

  return i;
}


//-------------------------------------------------------------------------------------------------------
// nextPow2
//-------------------------------------------------------------------------------------------------------
unsigned int nextPow2( unsigned int x ) {
  --x;
  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;
  return ++x;
}


//-------------------------------------------------------------------------------------------------------
// FreeDevice
//-------------------------------------------------------------------------------------------------------
void FreeDevice(void){

  cudaFree(d_Ra);
  cudaFree(d_Rb);
  cudaFree(d_Rc);
  cudaFree(d_Ba);
  cudaFree(d_Bb);
  cudaFree(d_Bc);

  cudaFree(d_Jpx);
  cudaFree(d_Jpy);
  cudaFree(d_Jpz);

  cudaFree(d_Energy);

}


//-------------------------------------------------------------------------------------------------------
// FreeHost
//-------------------------------------------------------------------------------------------------------
void FreeHost(void){

  Free(h_a);
  Free(h_b);
  Free(h_c);

  Free(h_Jpx);
  Free(h_Jpy);
  Free(h_Jpz);

  Free(h_Energy);

}


//-------------------------------------------------------------------------------------------------------
// Dump_Spin
//-------------------------------------------------------------------------------------------------------
void Dump_Spin(int count){

  FILE* tempfile;
  char namefile[MAXSTRLEN];

  memset(namefile,'\0',MAXSTRLEN); // NON CONTROLLO SE TORNA BENE O MENO

  if(count>=0){

    snprintf(namefile,sizeof(namefile),"spin_dump_%d",count);

    CUDA_SAFE_CALL( cudaMemcpy( h_a, d_Ra, ms_s, cudaMemcpyDeviceToHost ) );
    CUDA_SAFE_CALL( cudaMemcpy( h_b, d_Rb, ms_s, cudaMemcpyDeviceToHost ) );
    CUDA_SAFE_CALL( cudaMemcpy( h_c, d_Rc, ms_s, cudaMemcpyDeviceToHost ) );

  } else {
    snprintf(namefile,sizeof(namefile),"spin_dump_start");
  }

  tempfile = Fopen(namefile,"w");

  if(sizeof(REAL)==sizeof(float)) {

    for(int i = 0; i<ns_j; i++) fprintf(tempfile,"%f %f %f\n",h_a[i],h_b[i],h_c[i]);

  } else {

    for(int i = 0; i<ns_j; i++) fprintf(tempfile,"%lf %lf %lf\n",h_a[i],h_b[i],h_c[i]);

  }

  fflush(tempfile);

  fclose(tempfile);

}


//-------------------------------------------------------------------------------------------------------
// Dump_J
//-------------------------------------------------------------------------------------------------------
void Dump_J(void){

  FILE* tempfile;
  tempfile = Fopen("j_dump_start","w");

  if(sizeof(REAL)==sizeof(float)) {

    for(int i = 0; i<ns_j; i++) fprintf(tempfile,"%f %f %f\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);

  } else {

    for(int i = 0; i<ns_j; i++) fprintf(tempfile,"%lf %lf %lf\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);

  }


  fflush(tempfile);

  fclose(tempfile);

}


//-------------------------------------------------------------------------------------------------------
// Initialize_Spin
//-------------------------------------------------------------------------------------------------------
void Initialize_Spin(REAL* sa_a, REAL* sa_b, REAL* sa_c){

  REAL a, b, c;
  REAL invnorma, norma;

  for(int i = 0; i<h_nspin; i++){

    norma=two;

    while (norma>one){

      a=(RANDOM01)*two-one;
      b=(RANDOM01)*two-one;
      c=(RANDOM01)*two-one;

      norma=a*a+b*b+c*c;

    }

    invnorma=one/SQRT(norma);

    sa_a[i]=a*invnorma;
    sa_b[i]=b*invnorma;
    sa_c[i]=c*invnorma;

  }

}


//-------------------------------------------------------------------------------------------------------
// Read_Spin
//-------------------------------------------------------------------------------------------------------
void Read_Spin(char *inputfile, REAL* sr_a, REAL* sr_b, REAL* sr_c) {

  REAL a, b, c;
  FILE *fpi;
  fpi=Fopen(inputfile,"r");
  fprintf(stderr,"Reading spin from file %s\n",inputfile);

  for(int i = 0; i<ns_j; i++){
    if(sizeof(REAL)==sizeof(float)) {
      fscanf(fpi,"%f %f %f",&a,&b,&c);
    } else {
      fscanf(fpi,"%lf %lf %lf",&a,&b,&c);
    }
    sr_a[i]=a;
    sr_b[i]=b;
    sr_c[i]=c;
  }

  fclose(fpi);

}



//-------------------------------------------------------------------------------------------------------
// CopyDataOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyDataOnDevice(void){

  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Ra, ms_s ) );
  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Rb, ms_s ) );
  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Rc, ms_s ) );

  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Ba, ms_s ) );
  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Bb, ms_s ) );
  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Bc, ms_s ) );

  CUDA_SAFE_CALL( cudaMemcpy( d_Ra, h_a, ms_s, cudaMemcpyHostToDevice ) );
  CUDA_SAFE_CALL( cudaMemcpy( d_Rb, h_b, ms_s, cudaMemcpyHostToDevice ) );
  CUDA_SAFE_CALL( cudaMemcpy( d_Rc, h_c, ms_s, cudaMemcpyHostToDevice ) );

  CUDA_SAFE_CALL( cudaMemcpy( d_Ba, h_a, ms_s, cudaMemcpyHostToDevice ) );
  CUDA_SAFE_CALL( cudaMemcpy( d_Bb, h_b, ms_s, cudaMemcpyHostToDevice ) );
  CUDA_SAFE_CALL( cudaMemcpy( d_Bc, h_c, ms_s, cudaMemcpyHostToDevice ) );

  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Jpx, ms_j ) );
  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Jpy, ms_j ) );
  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Jpz, ms_j ) );

  CUDA_SAFE_CALL( cudaMemcpy( d_Jpx, h_Jpx, ms_j, cudaMemcpyHostToDevice ) );
  CUDA_SAFE_CALL( cudaMemcpy( d_Jpy, h_Jpy, ms_j, cudaMemcpyHostToDevice ) );
  CUDA_SAFE_CALL( cudaMemcpy( d_Jpz, h_Jpz, ms_j, cudaMemcpyHostToDevice ) );

  CUDA_SAFE_CALL( cudaMalloc( (void **) &d_Energy, ms_s ) );
  CUDA_SAFE_CALL( cudaMemset( d_Energy, 0, ms_s ) );
}


//-------------------------------------------------------------------------------------------------------
// Initialize_J
//-------------------------------------------------------------------------------------------------------
void Initialize_J(void){

  int j;

  switch (init_type_j){

  case 0:
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = ((RANDOM01)>zerofive) ? -1 : 1;
      h_Jpy[j] = ((RANDOM01)>zerofive) ? -1 : 1;
      h_Jpz[j] = ((RANDOM01)>zerofive) ? -1 : 1;
    }

    break;
  case 1:
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = h_Jpy[j] = h_Jpz[j] = one;
    }
    break;
  case -1:
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = h_Jpy[j] = h_Jpz[j] = -one;
    }

    break;
  case 2: //Caso Gaussiano
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = rand_gauss();
      h_Jpy[j] = rand_gauss();
      h_Jpz[j] = rand_gauss();
    }
    break;
  default: {
    writelog(TRUE,APPLICATION_RC,"Invalid initialization for J: %d", init_type_j);
    break;
  }
  }

}


//-------------------------------------------------------------------------------------------------------
// rand_gauss
//-------------------------------------------------------------------------------------------------------
REAL rand_gauss(void){

  REAL r,cosr;
  REAL a=one-RANDOM01;
  if(a<EPSILON) { a=EPSILON; }

  r=SQRT(-2.*LOG(a));

  cosr=COS((2*M_PI)*(RANDOM01));

  return r*cosr;

}


//-------------------------------------------------------------------------------------------------------
// Parameters_Reading
//-------------------------------------------------------------------------------------------------------
void Parameters_Reading(char *inputfile){

  dictionary *ini;
  char key[MAXSTRLEN];

  if(inputfile==NULL) { writelog(TRUE,APPLICATION_RC,"No inputfile"); }

  ini = iniparser_load(inputfile);

  if(ini==NULL) { writelog(TRUE,APPLICATION_RC,"Cannot parse file: %s\n", inputfile); }

  READINTFI(h_side,"Side");
  READINTFI(init_type_j,"Initialization type for J");
  READINTFI(nsteps,"Number of steps");
  READINTFI(overrelaxiter,"Number of overrelaxation iterations");
  READINTFI(overrelaxhits,"Number of overrelaxation hits");
  READINTFI(outfreq,"Output frequency");
  READINTFI(dumpfreq,"Dump configuration frequency");
  READREALFI(h_h,"Magnetic field");
  READREALFI(h_Beta,"Beta (1/T)");
  READINTFI(seed,"Random seed");

  DEVICEREADINTFI(numBlocks,"Number of blocks");
  DEVICEREADINTFI(blocksize_x,"Number of threads for block along X");
  DEVICEREADINTFI(numThreadsPerBlock,"Number of threads for block");

  DEVICEREADINTFI(RandNumBlocks,"RandGene num of blocks");
  DEVICEREADINTFI(RandNumThreadsPerBlock,"RandGene num of threads for block");

}


//-------------------------------------------------------------------------------------------------------
// Calculate_Sizes
//-------------------------------------------------------------------------------------------------------
void Calculate_Sizes(void){

  ns_j = V(h_side);

  ms_s = ns_j*sizeof(REAL);

  ms_j = ns_j*sizeof(REAL);

}


//-------------------------------------------------------------------------------------------------------
// getNumBlocksAndThreads
//-------------------------------------------------------------------------------------------------------
void getNumBlocksAndThreads(int n, int maxBlocks, int maxThreads, int &blocks,
                            int &threads) {
  if (n == 1)
    {
      threads = 1;
      blocks = 1;
    }
  else
    {
      threads = (n < maxThreads*2) ? nextPow2(n / 2) : maxThreads;
      blocks = max(1, n / (threads * 2));
    }

  blocks = min(maxBlocks, blocks);
}

//-------------------------------------------------------------------------------------------------------
// main
//-------------------------------------------------------------------------------------------------------
int main(int argc, char **argv){

  char *inputfile = NULL;
  char *po;
  int istep;
  int nblocks_red, nthreads_red;
  int cudadevice;
  struct cudaDeviceProp prop;
  int mpc, shmsize, mtpb, shmsidelen, shmitems;
  int number_of_blocks_x, number_of_blocks_y,
    number_of_blocks, number_of_threads_per_block;

  double ttimeover=0., ttimehb=0.;
  REAL normene;
  int j;
  char *readspin=NULL;

#if !defined(HOSTREDU)
  unsigned int retCnt;
#endif

  TIMER_DEF;

  for(int i = 1; i < argc; i++) {
    po = argv[i];
    if (*po++ == '-') {
      switch (*po++) {
      case '?':
        SKIPBLANK
          Usage(argv[0]);
        exit(OK);
        break;
      case 'v':
        SKIPBLANK
          verbose=TRUE;
        break;
      case 'i':
        SKIPBLANK
          inputfile=Strdup(po);
        break;
      case 'o':
        SKIPBLANK
          onlyOverrelaxation=true;
        break;
      case 'r':
        SKIPBLANK
          readspin=Strdup(po);
        break;
      case 'b':
        SKIPBLANK
          onlyHeatBath=true;
        break;

      default:
        Usage(argv[0]);
        exit(OK);
        break;
      }
    }
  }

  CUDA_SAFE_CALL(cudaGetDevice( &cudadevice ) );

  CUDA_SAFE_CALL(cudaGetDeviceProperties  (&prop, cudadevice) );
  mpc=prop.multiProcessorCount;
  mtpb=prop.maxThreadsPerBlock;
  shmsize=prop.sharedMemPerBlock;

  printf("Device %d: number of multiprocessors %d\n"
         "max number of threads per block %d; shared memory per block %d\n",
         cudadevice, mpc, mtpb, shmsize);

  Parameters_Reading(inputfile);

  shmsidelen=((shmsize-SHM_FOR_ARGUMENTS)/
              (NUMBER_OF_PLANES*NUMBER_OF_VARIABLES*sizeof(REAL))/
              (blocksize_x+TOTALHALO))-TOTALHALO;
  number_of_blocks_x=(h_side/blocksize_x);
  number_of_blocks_y=(h_side/shmsidelen)+(h_side%shmsidelen?1:0);

  number_of_blocks=number_of_blocks_x*number_of_blocks_y;

  if(number_of_blocks<mpc) {
    writelog(TRUE,APPLICATION_RC,"%d blocks are not enough! at least %d required\n",number_of_blocks, mpc);
  }

  number_of_threads_per_block=(blocksize_x*shmsidelen);
  printf("Number of blocks: %d\n",number_of_blocks);
  printf("Number of threads per block along x: %d\n",blocksize_x);
  printf("Number of threads per block along y: %d\n",shmsidelen);
  shmnitemspp=(shmsidelen+TOTALHALO)*(blocksize_x+TOTALHALO);
  printf("Area of shm plane: %d\n",shmnitemspp);
  shmitems=shmnitemspp*NUMBER_OF_PLANES*NUMBER_OF_VARIABLES;
  shmpervar=NUMBER_OF_PLANES*shmnitemspp;


  if (h_side<=0 || !pow_of_two(h_side)) {
    writelog(TRUE,APPLICATION_RC,"Invalid side: must be > 0 and power of two\n");
  }

  Calculate_Sizes();

  Init();
  if(readspin==NULL) {
    Initialize_Spin(h_a,h_b,h_c);
    Dump_Spin(-1);
  } else {
    Read_Spin(readspin,h_a,h_b,h_c);
  }

  Initialize_J();

  //      Dump_J();

  initRand(V(h_side));

  CopyConstOnDevice();

  CopyDataOnDevice();

  getNumBlocksAndThreads(h_nspin, numBlocks, numThreadsPerBlock,nblocks_red, nthreads_red);

  normene=one/(three*V(h_side));

  h_Energy[0]=zero;

  ComputeHEnergy(h_a, h_b, h_c, h_Jpx, h_Jpy, h_Jpz, h_Energy);

  printf("Initial Energy per spin on host: %f\n", normene*h_Energy[0]);

  CUDA_SAFE_CALL( cudaMalloc((void **)&r_odata,numBlocks*sizeof(REAL)) );

  dim3 dimGrid(numBlocks,1,1);
  dim3 dimBlock(numThreadsPerBlock,1,1);

  dim3 dimGridORHB(number_of_blocks_x, number_of_blocks_y,1);
  dim3 dimBlockORHB(blocksize_x, shmsidelen, 1);

  ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc, d_Jpx, d_Jpy, d_Jpz, d_Energy);

#if defined(HOSTREDU)

  CUDA_SAFE_CALL( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL),cudaMemcpyDeviceToHost) );

  h_totene=zeroEne;

  for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];

#else

  CUDA_SAFE_CALL( cudaThreadSynchronize() );
  retCnt=0;
  cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
  reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
  CUDA_SAFE_CALL( cudaThreadSynchronize() );
  CUDA_SAFE_CALL( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );

#endif

  printf("Initial Energy per spin: %f\n", normene*h_totene);

  int iterover=0;

  for(istep=0; istep<nsteps; istep++) {

    if(!onlyHeatBath){

      TIMER_START;

      for(j=0; j<overrelaxiter; j++, iterover++) {

        Overrelaxation<<< dimGridORHB,dimBlockORHB, shmitems*sizeof(REAL) >>>(d_Ra, d_Rb, d_Rc,
                                                                          d_Jpx, d_Jpy, d_Jpz,
                                                                          iterover&1);

        CUT_CHECK_ERROR("Overrelaxation() execution failed\n");

        CUDA_SAFE_CALL( cudaThreadSynchronize() );
      }

      TIMER_STOP;
      ttimeover+=TIMER_ELAPSED;

    }

    if(!onlyOverrelaxation) {

      TIMER_START;
      /* Here it goes the heat-bath */
      /* first of all, random number generation */
      /* then the heat-bath itself */

      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ONEONE,d_cRand);
      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_aRand);
      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_bRand);
      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROONE,d_Rand);

      if(istep==0&&FALSE){
        printRand(randNum,d_aRand,"daRand3");
        printRand(randNum,d_bRand,"dbRand3");
        printRand(randNum,d_cRand,"dcRand3");
        printRand(randNum,d_Rand,"dRand3");
        exit(1);
      }

      HeatBath<<< dimGridORHB,dimBlockORHB, shmitems*sizeof(REAL) >>>(d_Ra, d_Rb, d_Rc,
                                                                  d_Jpx, d_Jpy, d_Jpz,
                                                                  d_aRand,d_bRand,d_cRand,
                                                                  d_Rand,istep&1);

      CUT_CHECK_ERROR("HeatBath() execution failed\n");

      CUDA_SAFE_CALL( cudaThreadSynchronize() );


      TIMER_STOP;
      ttimehb+=TIMER_ELAPSED;

    }

    if(istep && (istep%outfreq)==0) {

      ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc,d_Jpx, d_Jpy, d_Jpz, d_Energy);

      CUDA_SAFE_CALL( cudaThreadSynchronize() );

#if defined(HOSTREDU)
      CUDA_SAFE_CALL( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL), cudaMemcpyDeviceToHost) );

      h_totene=zero;

      for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];

#else
      retCnt=0;
      cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
      reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
      CUDA_SAFE_CALL( cudaMemcpy(&h_totene, r_odata, sizeof(REAL), cudaMemcpyDeviceToHost) );
#endif
      printf("Energy per spin at step %d: %f\n", istep,normene*h_totene);

    }

    if(istep && ((istep%dumpfreq)==0)) {
      //Dump_Spin(istep);
    }

  }

  ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc, d_Jpx, d_Jpy, d_Jpz, d_Energy);

#if defined(HOSTREDU)

  CUDA_SAFE_CALL( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL), cudaMemcpyDeviceToHost) );

  h_totene=zero;

  for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];

#else
  CUDA_SAFE_CALL( cudaThreadSynchronize() );
  reduceSinglePass(ns_j, nthreads_red, nblocks_red, d_Energy, r_odata);
  CUDA_SAFE_CALL( cudaMemcpy(&h_totene, r_odata, sizeof(REAL), cudaMemcpyDeviceToHost) );
#endif

  printf("Final Energy per spin: %f\n",normene*h_totene);

  if(!onlyHeatBath){
    printf("Total time for %d iterations of overrelaxation on a %d^3 cube: %f seconds\n", nsteps*overrelaxiter*overrelaxhits, h_side, ttimeover*0.000001);

    printf("Time for spin: %f microseconds\n",ttimeover/(overrelaxiter*overrelaxhits*(double)nsteps*V(h_side)));
  }

  if(!onlyOverrelaxation){
    printf("Total time for %d iterations of heat bath on a %d^3 cube: %f seconds\n", nsteps, h_side, ttimehb*0.000001);

    printf("Time for spin: %f microseconds\n",ttimehb/((double)nsteps*V(h_side)));
  }


  Dump_Spin(nsteps);

  FreeDevice();

  FreeHost();

  FreeRand();

  return 0;

}


//-------------------------------------------------------------------------------------------------------
// Overrelaxation
//-------------------------------------------------------------------------------------------------------
#define MULT(a,b) (a)*(b) // __mul24((a),(b))

extern __shared__ REAL shmbuf[];
#define INDEX(i,j,j_off)  ((i) +MULT((j),(j_off)))
#define LOADPOINT(t,s)   sh_a[(t)]=a[(s)];      \
  sh_b[(t)]=b[(s)];                             \
  sh_c[(t)]=c[(s)];                             \
  sh_Jpx[(t)]=Jpx[(s)];                         \
  sh_Jpy[(t)]=Jpy[(s)];                         \
  sh_Jpz[(t)]=Jpz[(s)];


__global__ void Overrelaxation(REAL *a, REAL *b, REAL *c,
                               REAL *Jpx, REAL *Jpy, REAL *Jpz, int odditer) {

  unsigned int ix, iy, iz;
  REAL as, bs, cs;

  int spx, smx, spy, smy, spz, smz;

  REAL factor;

  int ig, ind, tbd, cpm1, cp, cpp1, cpp2, zoff, zoffp1, zoffp2, indxy;
  int halox=0, haloy=0, bulk=1, hshmix, hglmix, hshmiy, hglmiy, white=0;
  int hits;

  REAL *sh_a=shmbuf;
  REAL *sh_b=&sh_a[d_const_shmpervar];
  REAL *sh_c=&sh_b[d_const_shmpervar];
  REAL *sh_Jpx=&sh_c[d_const_shmpervar];
  REAL *sh_Jpy=&sh_Jpx[d_const_shmpervar];
  REAL *sh_Jpz=&sh_Jpy[d_const_shmpervar];

  ix=INDEX(threadIdx.x,blockIdx.x,blockDim.x);
  iy=INDEX(threadIdx.y,blockIdx.y,blockDim.y);

  tbd=(ix<d_side) && (iy<d_side);

#if defined(DEBUG)
  if(tbd) {
    printf("thread x: %d, thread y: %d, ix=%d, iy=%d, blockidX=%d, blockIdY=%d\n",threadIdx.x, threadIdx.y, ix, iy, blockIdx.x, blockIdx.y);
  }
#endif

  smz=ix+(iy<<d_sdbL)+(d_maskL<<d_sdbLs); //wrap around the lattice
  spz=ix+(iy<<d_sdbL)+(1<<d_sdbLs);

  cpm1=0;
  cp=1;
  cpp1=2;
  cpp2=3;

  ind=INDEX(threadIdx.x+1,threadIdx.y+1,blockDim.x+TOTALHALO);

  zoff=d_const_shmnitemspp;

  if(tbd) {
    if((iy&1)==(ix&1)) { white=1; }
    // load lower and current plane (with no halo)
    LOADPOINT(ind, smz);
    LOADPOINT(ind+zoff, ix+(iy<<d_sdbL));

    if(threadIdx.x==0) {
      hshmix=ind-1;
      hglmix=((ix-1)&d_maskL)+(iy<<d_sdbL);
      halox=1;
    }
    if(threadIdx.y==0) {
      hshmiy=ind-(blockDim.x+TOTALHALO);
      hglmiy=ix+(((iy-1)&d_maskL)<<d_sdbL);
      haloy=1;
    }
    if(threadIdx.x==(blockDim.x-1) || ix==(d_side-1)) {
      hshmix=ind+1;
      hglmix=((ix+1)&d_maskL)+(iy<<d_sdbL);
      halox=1;
    }
    if(threadIdx.y==(blockDim.y-1) || iy==(d_side-1)) {
      hshmiy=ind+(blockDim.x+TOTALHALO);
      hglmiy=ix+(((iy+1)&d_maskL)<<d_sdbL);
      haloy=1;
    }
    if(halox || haloy) {
      bulk=0;
    }

  }

  indxy=ind;

  for(iz=0; iz<d_side; iz+=2) {

    //lindex = ix + (iy * d_side) + (iz * d_side * d_side); //AGGIUNTO IO PER DEBUG

    //printf("%d\n",lindex); //AGGIUNTO IO PER DEBUG

    if(tbd) {

      zoff=MULT(cp,d_const_shmnitemspp);
      zoffp1=MULT(cpp1,d_const_shmnitemspp);
      zoffp2=MULT(cpp2,d_const_shmnitemspp);
      // load two planes
      ig=ix+(iy<<d_sdbL)+(((iz+1)&d_maskL)<<d_sdbLs);
      LOADPOINT(indxy+zoffp1, ig);

      ig=ix+(iy<<d_sdbL)+(((iz+2)&d_maskL)<<d_sdbLs);
      LOADPOINT(indxy+zoffp2, ig);

      // Load the halo of the two planes
      if(halox) {
        LOADPOINT(hshmix+zoff,hglmix+(iz<<d_sdbLs));
        LOADPOINT(hshmix+zoffp1,hglmix+(((iz+1)&d_maskL)<<d_sdbLs));
      }

      if(haloy) {
        LOADPOINT(hshmiy+zoff,hglmiy+(iz<<d_sdbLs));
        LOADPOINT(hshmiy+zoffp1,hglmiy+(((iz+1)&d_maskL)<<d_sdbLs));
      }

    }

    __syncthreads();
    for(hits=0; hits<d_overrelaxhits; hits++) {
      if(tbd) {

        ind=indxy+(white?zoff:zoffp1); //Add the offset along Z for the current plane
        spx=ind+1;
        as=sh_a[spx]*sh_Jpx[ind];
        bs=sh_b[spx]*sh_Jpx[ind];
        cs=sh_c[spx]*sh_Jpx[ind];

        smx=ind-1;
        as+=sh_a[smx]*sh_Jpx[smx];
        bs+=sh_b[smx]*sh_Jpx[smx];
        cs+=sh_c[smx]*sh_Jpx[smx];

        spy=ind+(blockDim.x+TOTALHALO);
        as+=sh_a[spy]*sh_Jpy[ind];
        bs+=sh_b[spy]*sh_Jpy[ind];
        cs+=sh_c[spy]*sh_Jpy[ind];

        smy=ind-(blockDim.x+TOTALHALO);
        as+=sh_a[smy]*sh_Jpy[smy];
        bs+=sh_b[smy]*sh_Jpy[smy];
        cs+=sh_c[smy]*sh_Jpy[smy];

        spz=ind+(white?(zoffp1-zoff):(zoffp2-zoffp1));
        as+=sh_a[spz]*sh_Jpz[ind];
        bs+=sh_b[spz]*sh_Jpz[ind];
        cs+=sh_c[spz]*sh_Jpz[ind];

        smz=ind+(white?(MULT(cpm1,d_const_shmnitemspp)-zoff):(zoff-zoffp1));
        as+=sh_a[smz]*sh_Jpz[smz];
        bs+=sh_b[smz]*sh_Jpz[smz];
        cs+=sh_c[smz]*sh_Jpz[smz];

        cs+=d_h;

        //Overrelaxation

        factor=d_const_two*(as*sh_a[ind]+bs*sh_b[ind]+cs*sh_c[ind])/ (as*as+bs*bs+cs*cs);

        if(bulk || odditer) {
#if defined(DEBUG)
          printf("%d %d %d %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f",ix,iy, iz, sh_a[ind],sh_b[ind],sh_c[ind], as,bs,cs, factor);
#endif
          sh_a[ind]=as*factor-sh_a[ind];
          sh_b[ind]=bs*factor-sh_b[ind];
          sh_c[ind]=cs*factor-sh_c[ind];
        }
      }

      __syncthreads();

      if(tbd) {

        ind=indxy+(white?zoffp1:zoff); //Add the offset along Z for the current plane
        spx=ind+1;
        as=sh_a[spx]*sh_Jpx[ind];
        bs=sh_b[spx]*sh_Jpx[ind];
        cs=sh_c[spx]*sh_Jpx[ind];

        smx=ind-1;
        as+=sh_a[smx]*sh_Jpx[smx];
        bs+=sh_b[smx]*sh_Jpx[smx];
        cs+=sh_c[smx]*sh_Jpx[smx];

        spy=ind+(blockDim.x+TOTALHALO);
        as+=sh_a[spy]*sh_Jpy[ind];
        bs+=sh_b[spy]*sh_Jpy[ind];
        cs+=sh_c[spy]*sh_Jpy[ind];

        smy=ind-(blockDim.x+TOTALHALO);
        as+=sh_a[smy]*sh_Jpy[smy];
        bs+=sh_b[smy]*sh_Jpy[smy];
        cs+=sh_c[smy]*sh_Jpy[smy];

        spz=ind+(white?(zoffp2-zoffp1):(zoffp1-zoff));
        as+=sh_a[spz]*sh_Jpz[ind];
        bs+=sh_b[spz]*sh_Jpz[ind];
        cs+=sh_c[spz]*sh_Jpz[ind];

        smz=ind+(white?(zoff-zoffp1):(MULT(cpm1,d_const_shmnitemspp)-zoff));
        as+=sh_a[smz]*sh_Jpz[smz];
        bs+=sh_b[smz]*sh_Jpz[smz];
        cs+=sh_c[smz]*sh_Jpz[smz];

        cs+=d_h;

        //Overrelaxation

        factor=d_const_two*(as*sh_a[ind]+bs*sh_b[ind]+cs*sh_c[ind])/(as*as+bs*bs+cs*cs);

        if(bulk || !odditer) {
#if defined(DEBUG)
          printf("%d %d %d %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f",ix,iy, iz, sh_a[ind],sh_b[ind],sh_c[ind], as,bs,cs, factor);
#endif
          sh_a[ind]=as*factor-sh_a[ind];
          sh_b[ind]=bs*factor-sh_b[ind];
          sh_c[ind]=cs*factor-sh_c[ind];
        }
      }
      __syncthreads();
    }
    if(tbd) {
      ig=ix+(iy<<d_sdbL)+(iz<<d_sdbLs);
      ind=indxy+zoff; //Add the offset along Z for the current plane

      a[ig]=sh_a[ind];
      b[ig]=sh_b[ind];
      c[ig]=sh_c[ind];

#if defined(DEBUG)
      printf("%d %7.5f %7.5f %7.5f\n", ig, a[ig],b[ig],c[ig]);
#endif

      ig=ix+(iy<<d_sdbL)+((iz+1)<<d_sdbLs);
      ind=indxy+zoffp1; //Add the offset along Z for the plus1 plane

      a[ig]=sh_a[ind];
      b[ig]=sh_b[ind];
      c[ig]=sh_c[ind];

#if defined(DEBUG)
      printf("%d %7.5f %7.5f %7.5f\n", ig, a[ig],b[ig],c[ig]);
#endif
    }

    cp+=2;   if(cp>=NUMBER_OF_PLANES)   {cp-=NUMBER_OF_PLANES;}
    cpm1+=2; if(cpm1>=NUMBER_OF_PLANES) {cpm1-=NUMBER_OF_PLANES;}
    cpp1+=2; if(cpp1>=NUMBER_OF_PLANES) {cpp1-=NUMBER_OF_PLANES;}
    cpp2+=2; if(cpp2>=NUMBER_OF_PLANES) {cpp2-=NUMBER_OF_PLANES;}

  }

}

//-------------------------------------------------------------------------------------------------------
// ComputeEnergy
//-------------------------------------------------------------------------------------------------------
__global__ void ComputeEnergy(REAL *a, REAL *b, REAL *c,
                              REAL *Jpx, REAL *Jpy, REAL *Jpz,
                              REAL *energy) {

  unsigned int ix, iy, iz;
  int spx, spy, spz;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x;


  for(i=istart+tid; i<d_nspin; i+=num_threads) {
    iz=i>>d_sdbLs;
    iy=(i - (iz<<d_sdbLs))>>d_sdbL;
    ix=(i - (iz<<d_sdbLs) - (iy<<d_sdbL));

    spx=((ix+1)&d_maskL)+(iy<<d_sdbL)+(iz<<d_sdbLs);
    spy=ix+(((iy+1)&d_maskL)<<d_sdbL)+(iz<<d_sdbLs);
    spz=ix+(iy<<d_sdbL)+(((iz+1)&d_maskL)<<d_sdbLs);

    energy[i]=a[i]*(a[spx]*Jpx[i]+
                    a[spy]*Jpy[i]+
                    a[spz]*Jpz[i])+
      b[i]*(b[spx]*Jpx[i]+
            b[spy]*Jpy[i]+
            b[spz]*Jpz[i])+
      c[i]*(c[spx]*Jpx[i]+
            c[spy]*Jpy[i]+
            c[spz]*Jpz[i]);

  }

}


//-------------------------------------------------------------------------------------------------------
// ComputeHEnergy
//-------------------------------------------------------------------------------------------------------
void ComputeHEnergy(REAL *a, REAL *b, REAL *c,
                    REAL *Jpx, REAL *Jpy, REAL *Jpz,
                    REAL *energy) {

  unsigned int ix, iy, iz;
  int spx, spy, spz;

  REAL e;

  int i;

  for(i=0; i<h_nspin; i++) {
    iz=i>>h_sdbLs;
    iy=(i - (iz<<h_sdbLs))>>h_sdbL;
    ix=(i - (iz<<h_sdbLs) - (iy<<h_sdbL));

    spx=((ix+1)&h_maskL)+(iy<<h_sdbL)+(iz<<h_sdbLs);
    spy=ix+(((iy+1)&h_maskL)<<h_sdbL)+(iz<<h_sdbLs);
    spz=ix+(iy<<h_sdbL)+(((iz+1)&h_maskL)<<h_sdbLs);

    e=a[i]*(a[spx]*Jpx[i]+
            a[spy]*Jpy[i]+
            a[spz]*Jpz[i])+
      b[i]*(b[spx]*Jpx[i]+
            b[spy]*Jpy[i]+
            b[spz]*Jpz[i])+
      c[i]*(c[spx]*Jpx[i]+
            c[spy]*Jpy[i]+
            c[spz]*Jpz[i]);
    energy[0]+=e;

  }

}


//-------------------------------------------------------------------------------------------------------
// ComputeRandomGPU
//-------------------------------------------------------------------------------------------------------
void ComputeRandomGPU(int numBlocksRand, int numThreadsPerBlockRand, int mode, REAL* array){

  seedRand(0);

  CUDA_SAFE_CALL( cudaThreadSynchronize() );

  dim3 dimGrid(numBlocksRand,1,1);

  dim3 dimBlock(numThreadsPerBlockRand,1,1);

  RandomGPU<<<dimGrid, numThreadsPerBlockRand>>>(array,mode);

  CUDA_SAFE_CALL( cudaThreadSynchronize() );

}


//-------------------------------------------------------------------------------------------------------
// FreeRand
//-------------------------------------------------------------------------------------------------------
void FreeRand(void){

  CUDA_SAFE_CALL( cudaFree(d_aRand) );
  CUDA_SAFE_CALL( cudaFree(d_bRand) );
  CUDA_SAFE_CALL( cudaFree(d_cRand) );
  CUDA_SAFE_CALL( cudaFree(d_Rand) );

}


//-------------------------------------------------------------------------------------------------------
// initRand
//-------------------------------------------------------------------------------------------------------
int initRand(int randomNumber){


  int n_per_rng;
  // int randNum; DA DISCOMMENTARE DOPO PROVA MESSO GLOBAL

  n_per_rng = iAlignUp(iDivUp(randomNumber, MT_RNG_COUNT), 2);

  randNum = MT_RNG_COUNT * n_per_rng;

  CUDA_SAFE_CALL( cudaMalloc((void **)&d_aRand, randNum * sizeof(REAL)) );
  CUDA_SAFE_CALL( cudaMalloc((void **)&d_bRand, randNum * sizeof(REAL)) );
  CUDA_SAFE_CALL( cudaMalloc((void **)&d_cRand, randNum * sizeof(REAL)) );
  CUDA_SAFE_CALL( cudaMalloc((void **)&d_Rand, randNum * sizeof(REAL)) );

  CUDA_SAFE_CALL( cudaMemcpyToSymbol("NPerRng",&n_per_rng,
                                     sizeof(int),0,cudaMemcpyHostToDevice) );

  return randNum;

}


//-------------------------------------------------------------------------------------------------------
// seedRand
//-------------------------------------------------------------------------------------------------------
void seedRand(int seed){

  if(seed==0) seed=(int)(drand48()*1023); // UNA SCHIFEZZA

  static mt_struct_stripped h_MT[MT_RNG_COUNT];
  static mt_struct MT[MT_RNG_COUNT];

  FILE *fd = Fopen("./data/MersenneTwister.raw", "rb");

  for (int i = 0; i < MT_RNG_COUNT; i++){
    //Inline structure size for compatibility,
    //since pointer types are 8-byte on 64-bit systems (unused *state variable)
    Fread(MT + i, 16 /* sizeof(mt_struct) */ * sizeof(int), 1, fd);

  }

  fclose(fd);

  fd = Fopen("./data/MersenneTwister.dat", "rb");

  Fread(h_MT, sizeof(h_MT), 1, fd);

  fclose(fd);

  //Need to be thread-safe
  mt_struct_stripped *MT_temp = (mt_struct_stripped *)Malloc(MT_RNG_COUNT * sizeof(mt_struct_stripped));

  for(int i = 0; i < MT_RNG_COUNT; i++){
    MT_temp[i] = h_MT[i];
    MT_temp[i].seed = seed;
  }

  CUDA_SAFE_CALL( cudaMemcpyToSymbol("ds_MT",MT_temp,sizeof(h_MT),0,cudaMemcpyHostToDevice) );

  free(MT_temp);

}


//-------------------------------------------------------------------------------------------------------
// printRand
//-------------------------------------------------------------------------------------------------------
void printRand(int randomNumber,REAL* array,char* namefile){

  REAL * h_Rand  = (REAL *)malloc(randomNumber * sizeof(REAL));

  CUDA_SAFE_CALL( cudaMemcpy(h_Rand, array, randomNumber * sizeof(REAL), cudaMemcpyDeviceToHost) );

  FILE* out = Fopen(namefile, "w");

  for(int x=0;x<randomNumber;x++) { fprintf(out,"%f\n",h_Rand[x]); }

  fclose(out);

}



////////////////////////////////////////////////////////////////////////////////
// Write MT_RNG_COUNT vertical lanes of NPerRng random numbers to *d_Random.
// For coalesced global writes MT_RNG_COUNT should be a multiple of warp size.
// Initial states for each generator are the same, since the states are
// initialized from the global seed. In order to improve distribution properties
// on small NPerRng supply dedicated (local) seed to each twister.
// The local seeds, in their turn, can be extracted from global seed
// by means of any simple random number generator, like LCG.
////////////////////////////////////////////////////////////////////////////////
__global__ void RandomGPU( REAL *d_Random, int mode ){

  const int tid = blockDim.x * blockIdx.x + threadIdx.x;
  const int THREAD_N = blockDim.x * gridDim.x;

  int iState, iState1, iStateM, iOut;
  unsigned int mti, mti1, mtiM, x;
  unsigned int mt[MT_NN];

  REAL temp_rand;

  for(int iRng = tid; iRng < MT_RNG_COUNT; iRng += THREAD_N){
    //Load bit-vector Mersenne Twister parameters
    mt_struct_stripped config = ds_MT[iRng];

    //Initialize current state
    mt[0] = config.seed;
    for(iState = 1; iState < MT_NN; iState++)
      mt[iState] = (1812433253U * (mt[iState - 1] ^ (mt[iState - 1] >> 30)) + iState) & MT_WMASK;

    iState = 0;
    mti1 = mt[0];
    for(iOut = 0; iOut < NPerRng; iOut++){
      //iState1 = (iState +     1) % MT_NN
      //iStateM = (iState + MT_MM) % MT_NN
      iState1 = iState + 1;
      iStateM = iState + MT_MM;
      if(iState1 >= MT_NN) iState1 -= MT_NN;
      if(iStateM >= MT_NN) iStateM -= MT_NN;
      mti  = mti1;
      mti1 = mt[iState1];
      mtiM = mt[iStateM];

      x    = (mti & MT_UMASK) | (mti1 & MT_LMASK);
      x    =  mtiM ^ (x >> 1) ^ ((x & 1) ? config.matrix_a : 0);
      mt[iState] = x;
      iState = iState1;

      //Tempering transformation
      x ^= (x >> MT_SHIFT0);
      x ^= (x << MT_SHIFTB) & config.mask_b;
      x ^= (x << MT_SHIFTC) & config.mask_c;
      x ^= (x >> MT_SHIFT1);

      temp_rand = ((REAL)x + 1.0f) / 4294967296.0f; //Convert to (0, 1] float


      if(mode==ZEROONE) { // write to global memory
        d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand;
      } else if(mode==ONEONE) {
        d_Random[iRng + iOut * MT_RNG_COUNT] = 2.0f*temp_rand-1.0f;
      }else {
        d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand - 0.5f;
      }
    }
  }
}



//-------------------------------------------------------------------------------------------------------
// HeatBath
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBath(REAL *a, REAL *b, REAL *c,
                         REAL *Jpx, REAL *Jpy, REAL *Jpz,
                         REAL *atry, REAL *btry, REAL *ctry,
                         REAL *random, int odditer) {

  unsigned int ix, iy, iz;

  int skip;

  REAL as, bs, cs;
  REAL anew, bnew, cnew;
  REAL hk, bhk, projection, pl, invnorm;

  int spx, smx, spy, smy, spz, smz;

  REAL invexpbhk2;

  REAL factor;

  int ig, ind, tbd, cpm1, cp, cpp1, cpp2, zoff, zoffp1, zoffp2, indxy;
  int halox=0, haloy=0, bulk=1, hshmix, hglmix, hshmiy, hglmiy, white=0;

  REAL *sh_a=shmbuf;
  REAL *sh_b=&sh_a[d_const_shmpervar];
  REAL *sh_c=&sh_b[d_const_shmpervar];
  REAL *sh_Jpx=&sh_c[d_const_shmpervar];
  REAL *sh_Jpy=&sh_Jpx[d_const_shmpervar];
  REAL *sh_Jpz=&sh_Jpy[d_const_shmpervar];

  ix=INDEX(threadIdx.x,blockIdx.x,blockDim.x);
  iy=INDEX(threadIdx.y,blockIdx.y,blockDim.y);

  tbd=(ix<d_side) && (iy<d_side);

  int lindex;

#if defined(DEBUG)
  if(tbd) {
    printf("thread x: %d, thread y: %d, ix=%d, iy=%d, blockidX=%d, blockIdY=%d\n",threadIdx.x, threadIdx.y, ix, iy, blockIdx.x, blockIdx.y);
  }
#endif

  smz=ix+(iy<<d_sdbL)+(d_maskL<<d_sdbLs); //wrap around the lattice
  spz=ix+(iy<<d_sdbL)+(1<<d_sdbLs);

  cpm1=0;
  cp=1;
  cpp1=2;
  cpp2=3;

  ind=INDEX(threadIdx.x+1,threadIdx.y+1,blockDim.x+TOTALHALO);

  zoff=d_const_shmnitemspp;

  if(tbd) {
    if((iy&1)==(ix&1)) { white=1; }
    // load lower and current plane (with no halo)
    LOADPOINT(ind, smz);
    LOADPOINT(ind+zoff, ix+(iy<<d_sdbL));

    if(threadIdx.x==0) {
      hshmix=ind-1;
      hglmix=((ix-1)&d_maskL)+(iy<<d_sdbL);
      halox=1;
    }
    if(threadIdx.y==0) {
      hshmiy=ind-(blockDim.x+TOTALHALO);
      hglmiy=ix+(((iy-1)&d_maskL)<<d_sdbL);
      haloy=1;
    }
    if(threadIdx.x==(blockDim.x-1) || ix==(d_side-1)) {
      hshmix=ind+1;
      hglmix=((ix+1)&d_maskL)+(iy<<d_sdbL);
      halox=1;
    }
    if(threadIdx.y==(blockDim.y-1) || iy==(d_side-1)) {
      hshmiy=ind+(blockDim.x+TOTALHALO);
      hglmiy=ix+(((iy+1)&d_maskL)<<d_sdbL);
      haloy=1;
    }
    if(halox || haloy) {
      bulk=0;
    }

  }

  indxy=ind;

  for(iz=0; iz<d_side; iz+=2) {
    if(tbd) {
      zoff=MULT(cp,d_const_shmnitemspp);
      zoffp1=MULT(cpp1,d_const_shmnitemspp);
      zoffp2=MULT(cpp2,d_const_shmnitemspp);
      // load two planes
      ig=ix+(iy<<d_sdbL)+(((iz+1)&d_maskL)<<d_sdbLs);
      LOADPOINT(indxy+zoffp1, ig);

      ig=ix+(iy<<d_sdbL)+(((iz+2)&d_maskL)<<d_sdbLs);
      LOADPOINT(indxy+zoffp2, ig);

      // Load the halo of the two planes
      if(halox) {
        LOADPOINT(hshmix+zoff,hglmix+(iz<<d_sdbLs));
        LOADPOINT(hshmix+zoffp1,hglmix+(((iz+1)&d_maskL)<<d_sdbLs));
      }
      if(haloy) {
        LOADPOINT(hshmiy+zoff,hglmiy+(iz<<d_sdbLs));
        LOADPOINT(hshmiy+zoffp1,hglmiy+(((iz+1)&d_maskL)<<d_sdbLs));
      }

      ind=indxy+(white?zoff:zoffp1); //Add the offset along Z for the current plane
      spx=ind+1;
      smx=ind-1;
      spy=ind+(blockDim.x+TOTALHALO);
      smy=ind-(blockDim.x+TOTALHALO);
      spz=ind+(white?(zoffp1-zoff):(zoffp2-zoffp1));
      smz=ind+(white?(MULT(cpm1,d_const_shmnitemspp)-zoff):(zoff-zoffp1));
    }
    __syncthreads();
    if(tbd) {
      lindex = ix + (iy * d_side) + (iz * d_side * d_side);

      anew=atry[lindex];
      bnew=btry[lindex];
      pl=anew*anew+bnew*bnew;

      if(pl>d_const_zerotwentyfive) skip=TRUE;
      else skip=FALSE;

      cnew=ctry[lindex];

      as=sh_a[spx]*sh_Jpx[ind]+
        sh_a[spy]*sh_Jpy[ind]+
        sh_a[spz]*sh_Jpz[ind]+
        sh_a[smx]*sh_Jpx[smx]+
        sh_a[smy]*sh_Jpy[smy]+
        sh_a[smz]*sh_Jpz[smz];

      bs=sh_b[spx]*sh_Jpx[ind]+
        sh_b[spy]*sh_Jpy[ind]+
        sh_b[spz]*sh_Jpz[ind]+
        sh_b[smx]*sh_Jpx[smx]+
        sh_b[smy]*sh_Jpy[smy]+
        sh_b[smz]*sh_Jpz[smz];

      cs=sh_c[spx]*sh_Jpx[ind]+
        sh_c[spy]*sh_Jpy[ind]+
        sh_c[spz]*sh_Jpz[ind]+
        sh_c[smx]*sh_Jpx[smx]+
        sh_c[smy]*sh_Jpy[smy]+
        sh_c[smz]*sh_Jpz[smz];

      cs+=d_h;

      hk=SQRT(as*as+bs*bs+cs*cs);
      bhk=d_Beta*hk;

      invexpbhk2=EXP(-d_const_two*bhk);
      projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[lindex]))/bhk;
      if(projection<-d_const_one) { projection=-d_const_one; }

      invnorm=SQRT((d_const_one-cnew*cnew)/pl);

      anew*=invnorm;
      bnew*=invnorm;

      factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);

      anew = anew-factor*as;
      bnew = bnew-factor*bs;
      cnew = cnew-factor*cs;

      factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
      projection/=hk;

      if((bulk || odditer) && !skip) {
#if defined(DEBUG)
        printf("%d %d %d %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f",ix,iy, iz, sh_a[ind],sh_b[ind],sh_c[ind], as,bs,cs, factor);
#endif
        sh_a[ind]=as*projection+anew*factor;
        sh_b[ind]=bs*projection+bnew*factor;
        sh_c[ind]=cs*projection+cnew*factor;
      }
    }
    __syncthreads();

    if(tbd) {

      lindex = ix + (iy * d_side) + ((iz+1) * d_side * d_side);

      anew=atry[lindex];
      bnew=btry[lindex];
      pl=anew*anew+bnew*bnew;

      if(pl>d_const_zerotwentyfive) skip=TRUE;
      else skip=FALSE;

      cnew=ctry[lindex];

      ind=indxy+(white?zoffp1:zoff); //Add the offset along Z for the current plane
      spx=ind+1;
      smx=ind-1;
      spy=ind+(blockDim.x+TOTALHALO);
      smy=ind-(blockDim.x+TOTALHALO);
      spz=ind+(white?(zoffp2-zoffp1):(zoffp1-zoff));
      smz=ind+(white?(zoff-zoffp1):(MULT(cpm1,d_const_shmnitemspp)-zoff));

      as=sh_a[spx]*sh_Jpx[ind]+
        sh_a[spy]*sh_Jpy[ind]+
        sh_a[spz]*sh_Jpz[ind]+
        sh_a[smx]*sh_Jpx[smx]+
        sh_a[smy]*sh_Jpy[smy]+
        sh_a[smz]*sh_Jpz[smz];

      bs=sh_b[spx]*sh_Jpx[ind]+
        sh_b[spy]*sh_Jpy[ind]+
        sh_b[spz]*sh_Jpz[ind]+
        sh_b[smx]*sh_Jpx[smx]+
        sh_b[smy]*sh_Jpy[smy]+
        sh_b[smz]*sh_Jpz[smz];

      cs=sh_c[spx]*sh_Jpx[ind]+
        sh_c[spy]*sh_Jpy[ind]+
        sh_c[spz]*sh_Jpz[ind]+
        sh_c[smx]*sh_Jpx[smx]+
        sh_c[smy]*sh_Jpy[smy]+
        sh_c[smz]*sh_Jpz[smz];

      cs+=d_h;

      hk=SQRT(as*as+bs*bs+cs*cs);
      bhk=d_Beta*hk;

      invexpbhk2=EXP(-d_const_two*bhk);
      projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[lindex]))/bhk;
      if(projection<-d_const_one) { projection=-d_const_one; }

      invnorm=SQRT((d_const_one-cnew*cnew)/pl);

      anew*=invnorm;
      bnew*=invnorm;

      factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);

      anew = anew-factor*as;
      bnew = bnew-factor*bs;
      cnew = cnew-factor*cs;

      factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
      projection/=hk;

      if((bulk || !odditer) && !skip) {
#if defined(DEBUG)
        printf("%d %d %d %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f",ix,iy, iz, sh_a[ind],sh_b[ind],sh_c[ind], as,bs,cs, factor);
#endif
        sh_a[ind]=as*projection+anew*factor;
        sh_b[ind]=bs*projection+bnew*factor;
        sh_c[ind]=cs*projection+cnew*factor;
      }
      ig=ix+(iy<<d_sdbL)+(iz<<d_sdbLs);
      ind=indxy+zoff; //Add the offset along Z for the current plane

      a[ig]=sh_a[ind];
      b[ig]=sh_b[ind];
      c[ig]=sh_c[ind];

#if defined(DEBUG)
      printf("%d %7.5f %7.5f %7.5f\n", ig, a[ig],b[ig],c[ig]);
#endif

      ig=ix+(iy<<d_sdbL)+((iz+1)<<d_sdbLs);
      ind=indxy+zoffp1; //Add the offset along Z for the plus1 plane

      a[ig]=sh_a[ind];
      b[ig]=sh_b[ind];
      c[ig]=sh_c[ind];

#if defined(DEBUG)
      printf("%d %7.5f %7.5f %7.5f\n", ig, a[ig],b[ig],c[ig]);
#endif
    }

    cp+=2;   if(cp>=NUMBER_OF_PLANES)   {cp-=NUMBER_OF_PLANES;}
    cpm1+=2; if(cpm1>=NUMBER_OF_PLANES) {cpm1-=NUMBER_OF_PLANES;}
    cpp1+=2; if(cpp1>=NUMBER_OF_PLANES) {cpp1-=NUMBER_OF_PLANES;}
    cpp2+=2; if(cpp2>=NUMBER_OF_PLANES) {cpp2-=NUMBER_OF_PLANES;}

  }

}


